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Modelling  sheet-flow  sediment  transport  in 
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Sediment  transport  in  oscillatory  boundary  layers  is  a  process  that  drives  coastal 
geomorphological  change.  Most  formulae  for  bed-load  transport  in  nearshore  regions 
subsume  the  smallest-scale  physics  of  the  phenomena  by  parametrizing  interac¬ 
tions  amongst  particles.  In  contrast,  we  directly  simulate  granular  physics  in  the 
wave-bottom  boundary  layer  using  a  discrete-element  model  comprised  of  a  three- 
dimensional  particle  phase  coupled  to  a  one-dimensional  fluid  phase  via  Newton’s 
third  law  through  forces  of  buoyancy,  drag  and  added  mass.  The  particulate  sedi¬ 
ment  phase  is  modelled  using  discrete  particles  formed  to  approximate  natural  grains 
by  overlapping  two  spheres.  Both  the  size  of  each  sphere  and  the  degree  of  overlap  can 
be  varied  for  these  composite  particles  to  generate  a  range  of  non-spherical  grains. 
Simulations  of  particles  having  a  range  of  shapes  showed  that  the  critical  angle — the 
angle  at  which  a  grain  pile  will  fail  when  tilted  slowly  from  rest — increases  from 
approximately  26°  for  spherical  particles  to  nearly  39°  for  highly  non-spherical  com¬ 
posite  particles  having  a  dumbbell  shape.  Simulations  of  oscillatory  sheet  flow  were 
conducted  using  composite  particles  with  an  angle  of  repose  of  approximately  33° 
and  a  Corey  shape  factor  greater  than  about  0.8,  similar  to  the  properties  of  beach 
sand.  The  results  from  the  sheet-flow  simulations  with  composite  particles  agreed 
more  closely  with  laboratory  measurements  than  similar  simulations  conducted  using 
spherical  particles.  The  findings  suggest  that  particle  shape  may  be  an  important 
factor  for  determining  bed-load  flux,  particularly  for  larger  bed  slopes. 

Keywords:  sediment  transport;  bed  load;  wave-bottom  boundary  layer; 
discrete-element  model;  non-spherical  particle 


1.  Introduction 

Predicting  the  morphological  evolution  of  sandy  coastlines  from  the  highest  uprush 
of  swash  offshore  to  the  edge  of  the  surf  zone  where  waves  begin  to  break  is  a  problem 

One  contribution  of  12  to  a  Theme  ‘Discrete-element  modelling:  methods  and  applications  in  the  environ¬ 
mental  sciences’. 
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with  social,  economic  and  scientific  significance.  Despite  the  apparent  accessibility  of 
the  phenomena  of  interest — namely,  the  motion  of  sand  under  the  forcing  of  waves 
and  currents — the  predictive  capability  of  existing  models  for  nearshore  evolution  is 
limited  to  qualitative  estimates  of  bulk  behaviours.  Given  adequate  forcing,  sediment 
transport  in  the  nearshore  responds  in  two  dominant  modes:  bed-load  transport, 
in  which  grains  collide,  slide,  bounce  and  roll  in  close  proximity  to  the  bed;  and 
suspended  load,  in  which  grains  are  lifted  from  the  bed  by  fluid  turbulence  and 
suspended  within  the  water  column.  We  focus  on  a  special  case  of  bed-load  transport 
called  ‘sheet  flow’,  which  is  characterized  by  high-concentration  flows  resulting  from 
intense  wave/current  forcing  where  the  local  bed  becomes  nominally  planar  and 
sediment  moves  in  a  thin  sheet  that  may  be  several  centimetres  thick  with  a  distinct 
upper  surface.  When  the  waves  and  currents  are  large  we  believe  sheet  flow  is  a 
primary  agent  of  bathymetric  evolution. 

Direct  measurements  of  bed-load  transport  rates  in  the  surf  zone  do  not  exist 
due  to  a  fundamental  lack  of  technology.  The  use  of  sediment  traps  and  bathymet¬ 
ric  surveys  are  two  of  the  best  methods  available  to  infer  bed-load  transport  rates. 
Also,  detailed  measurements  of  bulk  properties  in  the  sheet  layer  (e.g.  concentra¬ 
tion,  velocity,  fluctuation  energy,  stress)  are  nearly  impossible  to  obtain  with  present 
technologies. 

The  relatively  simple  geometry  of  sheet-flow  transport  makes  it  a  well-suited  prob¬ 
lem  for  study  with  a  discrete-element  model  (DEM).  Drake  &  Calantoni  (2001)  used 
a  DEM  to  study  sheet-flow  transport  in  the  surf  zone  where  sediment  particles  were 
modelled  as  spheres.  Although  the  Drake  &  Calantoni  (2001)  DEM  is  a  useful  tool 
for  studying  the  general  behaviour  of  sheet  flow,  recent  advances  in  computing  tech- 
nology  now  make  it  possible  to  simulate  non-spherical  particle  shapes  and  examine 
the  effect  of  shape  on  bed-load  transport  phenomena.  Consider  the  critical  angle,  </>, 
as  that  at  which  a  pile  of  grains  will  experience  failure  solely  under  the  influence 
of  gravity.  The  critical  angle  of  beach  sand  is  much  larger  than  that  of  frictional, 
spherical  particles  having  the  same  distribution  of  sizes  as  the  natural  sand.  We 
hypothesize  that  particle  shape  significantly  influences  the  static  bed  strength  of 
non-cohesive  sediments,  a  key  determinant  of  bed-load  flux,  and  test  this  hypothesis 
by  approximating  the  shape  of  sand  grains  with  non-spherical  composite  particles 
constructed  by  joining  two  partial  spheres  of  different  radii  (figure  1). 

Simulations  demonstrated  that  there  exists  a  composite  particle  shape  where  a 
bed  of  these  particles  exhibits  a  coefficient  of  static  friction  similar  to  beach  sand, 
tan  4>  «  0.63.  The  same  composite  particle  was  then  used  to  simulate  sheet-flow  trans¬ 
port  conditions  based  on  the  laboratory  data  of  King  (1991).  The  sheet-flow  simu¬ 
lations  showed  that  the  bed-load  flux  agreed  well  with  the  laboratory  measurements 
for  gently  sloping  beds,  |tan/?|  <  0.1.  At  larger  bed  slopes  and  higher  transport  rates 
the  simulations  suggested  that  particle  shape  could  be  an  important  parameter  for 
predicting  bed-load  flux. 


2.  DEM  for  sheet-flow  transport 

The  new  DEM  described  here  alters  the  model  of  Drake  &  Calantoni  (2001)  by 
allowing  the  shape  of  sediment  grains  to  be  simulated  with  non-spherical  compos¬ 
ite  particles.  The  DEM  comprises  a  three-dimensional  Lagrangian  particle  model 
coupled  to  a  one-dimensional  Eulerian  fluid  model.  The  particle  and  fluid  phases  are 
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Figure  1.  An  image  of  a  composite  particle  is  shown  where  the  finite  plane  represents  the  plane 
of  separation  of  the  two  partial  spheres.  The  initial  body  coordinate  system  chosen  to  solve  the 
inertia  tensor  problem  put  the  axis  of  symmetry  of  the  composite  particle  along  the  positive 
z-axis.  The  symmetry  of  the  particle  allows  for  the  orientation  of  the  x-  and  y- axes  to  be 
arbitrarily  specified.  Here  only  the  location  of  the  (x,  ?/)-plane  is  shown. 

coupled  at  every  model  time-step  via  Newton’s  third  law  through  fluid— particle  forces 
of  buoyancy,  drag  and  added  mass.  An  eddy- viscosity-based  model  determined  from  a 
mixing-length-governed  vertical  momentum  transfer  between  discrete  fluid  elements 
completes  a  self-consistent  physical  description  of  the  sheet-flow  process.  The  soft 
sphere  interactions  between  individual  particles  are  modelled  exactly  as  in  Drake  & 
Calantoni  (2001).  However,  the  shape  of  simulated  sand  grains  is  no  longer  restricted 
to  be  spherical. 


(a)  The  composite  particle 

The  motivation  for  choosing  to  construct  non-spherical  particles  in  this  manner 
was  largely  guided  by  the  fact  that  for  the  composite  particle  the  moments  of  inertia 
are  solvable  analytically  in  the  appropriate  coordinate  system.  Additionally,  efficient 
algorithms  for  detecting  contacts  between  individual  spheres  could  be  used  to  detect 
contacts  between  composite  particles.  The  following  subsections  provide  the  neces¬ 
sary  analytical  solutions  and  explain  how  the  simulations  compute  motions  for  the 
composite  particle. 

(i)  Volume  calculation 

A  cylindrical  coordinate  system  with  the  positive  z-axis  pointing  from  the  centre 
of  sphere  1  to  the  centre  of  sphere  2,  where  R\  >  i?2,  was  a  suitable  choice  of 
body  coordinates  for  performing  the  volume  integration  of  the  composite  particle 
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Figure  2.  A  composite  particle  in  plan  view  (maximum  projected  area).  The  origin  for  the 
coordinate  system  was  chosen  at  the  centre  of  the  disc  of  intersection  of  the  two  partial  spheres. 
The  symmetry  axis  of  the  composite  particle  is  chosen  to  coincide  with  the  2-axis  allowing  the 
two  partial  spheres  to  be  separated  by  the  plane  2  =  0.  Partial  sphere  1  with  radius  is  placed 
below  the  plane  of  separation,  2  =  0,  while  partial  sphere  2  with  radius  R2  is  placed  above  the 
plane  of  separation  with  the  positive  2-axis  pointing  from  the  centre  of  partial  sphere  1  to  the 
centre  of  partial  sphere  2.  The  angles  71  and  72  are  measured  down  from  the  positive  2-axis 
about  the  centre  of  each  respective  partial  sphere. 

(figure  2).  The  plane  z  =  0  separates  the  composite  particle  into  two  partial  spheres 
allowing  the  volume  integration  to  be  separated  into  a  sum  over  the  two  partial 
spheres,  V  =  V1+V2.  Using  the  cylindrical  body  coordinate  system  for  the  composite 
particle  defined  in  figure  2,  where  the  origin  is  placed  at  the  centre  of  the  disk  of 
intersection  of  the  two  partial  spheres,  the  volume  integral  becomes 

ry/ cos7i)2  /* 2 n  r 0 

V  =  /  rdr  I  dO  dz 

•'O  Jo  */— /?i  (I+COS71) 

ry/RZ-(z+R2  COS72)2  r2ir  pR2{  1 -cos 72) 

+l  Tdrl  d70  d2  <21) 

Performing  the  integration  yields 

V  =  37r[fl?(2  +  3cos7i  -  cos3 71)  +  R\{ 2  -  3cos72  +  cos3 72)],  (2.2) 

the  volume  of  the  composite  particle. 
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(ii)  Quantifying  the  shape  of  natural  particles 

Researchers  have  constructed  a  variety  of  methods  to  describe  the  shape  of  a  nat¬ 
ural  clastic  sedimentary  particle  simply  with  a  single  parameter  in  addition  to  some 
measure  of  the  particle  diameter.  Classifying  and  describing  the  shape  of  natural 
particles  has  often  been  motivated  by  the  desire  to  find  a  functional  relationship 
between  shape  and  the  settling  velocity  (Dietrich  1982).  Our  motivation,  conversely, 
was  to  construct  a  composite  particle  that  exhibited  the  same  critical  angle  (or  inter¬ 
nal  angle  of  static  friction)  as  nat rural  beach  sand.  The  shape  of  composite  particles 
we  constructed  was  quantified  with  the  Corey  shape  factor  (CSF)  (Corey  1949). 
The  CSF  is  not  the  only  way  to  quantify  the  shape  of  natural  particles:  the  Wadell 
sphericity  (Wadell  1932)  and  the  maximum  projection  sphericity  (MPS)  are  other 
popular  factors  for  sediment  shape  classification  (e.g.  Middleton  &  Southard  1984). 
The  physical  limits  of  the  composite  particle  ranged  from  a  sphere  to  a  perfect  dumb¬ 
bell  (R2/R1  =  1,  71  =  0,  72  =  tt).  For  this  range  of  shapes  the  CSF  assumes  values 
from  0.71  to  1,  the  Wadell  sphericity  assumes  values  from  0.63  to  1,  and  the  MPS 
assumes  values  from  0.79  to  1.  In  the  limit  where  R2  and  71  go  to  zero  the  composite 
particle  becomes  a  sphere  and  the  model  presented  here  is  identical  to  the  DEM 
of  Drake  &  Calantoni  (2001).  In  the  following  sections  we  describe  simulations  to 
determine  a  composite  particle  shape  having  the  same  critical  angle  as  beach  sand, 
and  a  second  set  of  simulations  of  oscillatory  sheet  flow  using  that  composite  shape 
for  inter-comparison  between  laboratory  experiments. 

(iii)  Centre  of  mass 

The  position  of  the  centre  of  mass  (CM)  for  each  composite  particle  lies  on  the  z- 
axis  because  of  axial  symmetry,  and  the  ^-coordinate  of  the  CM  of  each  partial  sphere 
was  found  by  integrating  to  obtain  (e.g.  Marion  &  Thornton  1988,  example  8.1) 

Zi  -  (3  +  1 cos  71  +  5  cos2  71  —  -jtj  cos4  71),  (2.3  a) 

Z2  =  7r^~  ( j  -  §  COS 72  +  5  COS2  72  -  Y5  cos4  72),  (2-3 b) 

where  Mi  and  M2  are  the  masses  of  partial  sphere  1  and  2,  respectively.  Now  the 
z-coordinate  of  the  CM  for  the  composite  particle  is  given  by 

-  Mrnwg(Mi2‘ + Mm)-  <2-4) 

(iv)  Principal  moments  of  inertia 

The  inertia  tensor  for  the  composite  particle  is  computed  using  the  body  coordinate 
system  with  the  integration  limits  on  the  volume  defined  above.  The  inertia  tensor  is 
calculated  considering  the  density  of  each  partial  sphere  to  be  constant  and  uniform. 
Performing  the  integrations  (e.g.  Goldstein  1980,  eqn  (5.8))  gives  a  diagonal  inertia 
tensor  of  the  form 

h  0  0\ 

/  =  o  h  0  ,  (2.5) 

\0  0  h) 
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where 

h  —  h 


and 

h  =  ^  1  [(1  +  cos7i)(8  +  7  COS71  -  7  cos27i  -  3cos37x  +  3  cos4  71)] 

+  2[(1  “  cos  72) (8  -  7  cos  72  —  7  cos2  72  +  3  cos3  72  +  3  cos4  72)]. 

(2.66) 


irpiR\ 


0Q  [(1  +  cos7i)(16  +  29cos7i  +  11  cos2  71  —  cos3  71  +  cos4  71)] 

+  -P62q  2 [(1  -  cos72)(16  -  29 cos 72  +  llcos2  72  +  cos3 72  +  cos4 72)] 

(2.6  a) 


(v)  Principal  moments  of  inertia  about  the  CM 

The  principal  moments  of  inertia  about  the  CM  are  required  to  solve  Euler’s 
equations  for  the  rotational  motion  of  a  rigid  body.  In  general,  the  CM  of  a  composite 
particle  does  not  coincide  with  the  origin  of  the  body  coordinate  system.  Applying 
the  parallel  axis  theorem,  the  principal  moments  of  inertia  about  the  CM  of  the 
composite  particle  shown  in  figure  2  become 

Ji  =  h  -  {Mi  +  M2)zcM,  (2.7 a) 

1*2  “  h  -  {M\  +  M2)zqM ,  (2.7  b) 

^3  =  /3-  (2.7c) 

(vi)  Contact  detection  and  contact  forces 

Use  of  spheres  to  construct  composite  particles  allows  use  of  highly  optimized  rou¬ 
tines  to  detect  contacts  between  spheres.  The  contact  detection  algorithm  employed 
in  this  model  is  based  on  the  method  of  Munjiza  &  Andrews  (1998).  The  normal  and 
tangential  forces  generated  at  the  contact  points  between  partial  sphere  pairs  are 
calculated  using  the  same  interaction  models  as  Drake  &  Calantoni  (2001).  Forces 
are  then  re-expressed  in  terms  of  the  normal  and  tangential  components  relative  to 
the  CM  of  the  respective  composite  particle.  For  each  individual  force  determined  to 
be  tangential  to  the  CM  of  a  composite  particle,  a  moment  arm  is  constructed  and 
the  torque  is  calculated.  The  sum  of  the  torque  on  each  composite  particle  is  then 
calculated  along  with  the  sum  of  the  force  through  the  CM. 


(vii)  Fluid-fluid  interactions 

The  fluid  phase  of  the  model  used  the  same  interactions  as  Drake  &  Calantoni 
(2001).  A  one-dimensional  eddy-viscosity  model  determined  from  a  mixing  length 
was  driven  by  a  time-varying  horizontal  pressure  gradient  with  sinusoidal  forcing, 

Foccos(utf),  (2.8) 
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where  uj  is  the  angular  frequency.  The  applied  horizontal  pressure  gradient  was  con¬ 
stant  across  the  entire  domain  and  acted  on  both  the  fluid  and  the  particles  embedded 
in  the  wave-bottom  boundary  layer  (WBBL).  Computationally,  the  fluid  model  was 
implemented  with  a  stack  of  discrete  fluid  slabs,  such  that  if  the  CM  of  a  parti¬ 
cle  was  located  within  a  slab,  then  all  fluid-particle  interactions  occurred  with  that 
slab.  Likewise,  the  total  volume  of  the  particles  located  in  a  fluid  slab  appropriately 
reduced  the  volume  of  fluid  contained  in  the  slab. 

(viii)  Fluid-particle  interactions 

Fluid-particle  forces  are  determined  by  considering  the  forces  on  a  spherical  par¬ 
ticle  with  the  equivalent  mass.  The  governing  equation  for  translational  motion  of  a 
composite  particle  is  given  by  (e.g.  Madsen  1991) 


+  \pC^A\u  -  us\(u  -  ms)  +  F«p,  (2.9) 

where  ps  and  p  are  the  particle  and  fluid  densities,  respectively,  Vs  is  the  particle 
volume,  g  is  the  acceleration  due  to  gravity,  us  and  u  are  the  particle  and  fluid  veloc¬ 
ities,  respectively,  and  A  is  the  projected  area  of  the  equivalent  spherical  particle. 
All  derivatives  are  evaluated  at  the  CM  of  the  particle  unless  specifically  noted.  The 
first  term  on  the  right-hand  side  represents  the  particle  buoyancy,  the  second  term  is 
the  horizontal  pressure  gradient  acting  on  the  particle,  and  the  third  term  represents 
the  added-mass  effect  with  the  coefficient  of  added  mass  cm  =  0.5  (Batchelor  1967). 
The  fourth  term  represents  the  particle  drag  force  with  the  drag  coefficient  C £> ,  given 
by  an  approximate  fit  to  the  empirical  drag  law  for  spheres,  including  a  correction 
c*,  based  on  local  particle  concentration  (e.g.  Richardson  &  Zaki  1954), 

C£  =  c*(24iie“1  +  4ife71/2  +  0.4),  (2.10) 

where 

c*  =  (l-c-ic2)-5/2,  (2.11) 

Res  is  the  relative  particle  Reynolds  number,  and  c  is  the  local  particle  concentration. 
Fluid-particle  forces  resulting  from  the  rotation  of  particles  are  ignored  along  with 
lift  forces.  The  final  term  in  equation  (2.9),  F#,  represents  the  sum  of  interparticle  . 
forces  acting  through  the  CM  of  the  particle. 

( b )  Simulations  with  composite  particles 
(i)  Critical- angle  simulations 

A  suite  of  simulations  was  performed  using  composite  particles  that  nearly  span  the 
allowable  range  of  the  CSF.  The  critical  angle,  </>,  was  measured  for  eight  different 
beds  of  250  identical  composite  particles,  formed  by  fixing  the  value  of  the  ratio 
of  R2/R1  and  then  allowing  the  value  of  72  (defined  in  figure  2)  to  vary  from  0 
to  7r  in  increments  of  it/ 8.  Here  the  critical  angle  was  taken  to  be  the  angle  at 
which  a  pile  of  grains  fails  under  the  force  of  gravity,  which  is  also  referred  to  as  the 
internal  angle  of  static  friction.  For  the  case  of  72  =  0,  the  composite  particle  reduces 
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to  a  sphere  with  radius  given  by  Rx.  Values  of  Rx  =  2.5  mm  and  R2  =  1.25  mm 
were  chosen  so  that  R2/R1  =  5.  Simulations  were  performed  for  particles  under  the 
influence  of  gravity  in  a  vacuum;  no  interstitial  fluid  was  modelled  in  the  critical-angle 
simulations.  Simulations  began  by  settling  the  particles  from  a  regular  lattice  onto  a 
smooth  basal  plane  having  a  row  of  fixed  spheres  attached  at  one  end  to  prevent  the 
wholesale  sliding  of  the  granular  assemblage.  Periodic  boundary  conditions  were  used 
in  the  two  lateral  directions.  The  bed  was  then  ‘tilted’  by  rotating  the  gravity  force 
vector  in  discrete  increments  at  a  rate  of  0.1°  s-1.  Once  failure  occurred,  the  gravity 
vector  was  rotated  back  in  the  same  manner  until  it  was  again  orthogonal  to  the 
bed.  Ancillary  tests  showed  that  simulating  one  second  of  real  time  after  changing 
the  direction  of  the  gravity  vector  was  sufficient  to  measure  the  critical  angle  with  a 
precision  of  0.1°. 

(ii)  Sheet-flow  simulations  of  King ’s  experiments 

The  experiments  of  King  (1991)  served  as  a  useful  benchmark  to  validate  the  model 
presented  here  for  sheet-flow  transport  of  coarse  sediments  under  waves.  King’s  oscil¬ 
latory  flow  tunnel  experiments  using  a  distribution  of  quartz  sediment  with  a  mean 
diameter  of  1.1  mm  provide  sediment  transport  rates  for  half-cycle  oscillatory  flows  in 
a  number  of  horizontal  and  tilted-bed  configurations.  The  sheet-flow  simulations  were 
performed  with  the  composite  particle  shape  that  demonstrated  a  critical  angle  simi¬ 
lar  to  beach  sand,  <f>  ~  33°  (Bagnold  1956).  Sheet-flow  simulations  with  spheres  were 
also  performed  and  were  shown  to  emphasize  the  significant  effect  of  particle  shape 
on  bulk  sediment  transport  rates  predicted  by  the  model.  The  set  of  experiments 
simulated  here  use  a  half-period  sine  wave  with  a  period  T  =  4.64  s,  and  a  maximum 
velocity  amplitude  of  1.18  ms-1,  over  a  range  of  bed  slopes,  -0.096  ^  tan f)  ^  0.16, 
where  0  is  the  angle  of  inclination  of  the  bed  in  the  direction  of  flow. 

The  simulations  used  1600  particles  where  every  composite  particle  had  the  same 
shape  with  a  distribution  of  sizes  similar  to  the  distribution  used  in  the  King  (1991) 
experiments.  The  mass  of  each  of  the  1600  unique  composite  particles  was  held 
constant  as  the  shape  was  changed  to  a  sphere,  ensuring  that  the  only  difference 
in  the  sheet-flow  simulations  with  spheres  and  composite  particles  was  the  particle 
shape.  The  density  of  particles  in  all  simulations  was  constant  and  equivalent  to  the 
density  of  quartz  sediment,  ps  =  2650  kg  m"3.  For  these  sheet-flow  simulations,  the 
physical  dimensions  of  the  computational  domain  were  approximately  2  cm  along 
the  flow  direction,  1  cm  transverse  to  the  flow,  and  4  cm  in  the  vertical.  A  stack  of 
40  fluid  slabs,  each  1  mm  thick,  was  used  to  model  the  fluid. 


3.  Simulation  results 

(a)  Critical- angle  simulations 

The  simulations  performed  demonstrated  that  there  exists  a  composite  particle  shape 
that  does  exhibit  a  critical  angle  similar  to  beach  sand.  The  critical  angle,  <j>,  is  plotted 
as  a  function  of  the  CSF  in  figure  3.  Each  point  is  the  average  value  of  10  simulations 
with  one  standard  deviation  plotted  above  and  below.  A  bed  of  identical  composite 
particles  with  the  ratio  R2/R\  =  \  and  72  =  5tt/8  having  CSF  =  0.84  exhibited  an 
average  critical  angle  of  32.9°.  Naturally  occurring  sediments  have  an  upper  mean 


Phil.  Trans.  R.  Soc.  Loud.  A  (2004) 


Sheet- flow  sediment  transport  in  wave-bottom  boundary  layers  04TA2210/9 


Figure  3.  Each  triangle  is  the  average  critical  angle  of  10  successive  simulations  with  one  standard 
deviation  plotted  above  and  below  each  point  for  the  ratio  R2/R1  =  with  the  value  of  72 
increasing  smoothly  from  0  to  7r,  in  increments  of  7t/8,  as  the  CSF  decreases  from  1  to  0.76. 
Shown  in  the  upper  panel  are  images  of  composite  particles  located  in  the  position  of  their 
corresponding  CSF  value  with  72  —  7r,  37t/4,  7t/4  and  tt/S  from  left  to  right.  The  circle  represents 
the  critical  angle  for  a  bed  of  identical  spheres,  CSF  =  1.  The  composite  particle  with  the 
ratio,  R2 / R\  —  72  =  57r/8,  having  CSF  =  0.84,  exhibited  an  average  critical  angle  of  32.9°, 

suggesting  that  for  the  range  of  shapes  simulated  this  composite  particle  possesses  properties 
very  similar  to  beach  sand. 

CSF  value  of  about  0.8  (Dietrich  1982),  suggesting  that  we  found  a  composite  particle 
that  simultaneously  satisfied  shape  and  critical-angle  characteristics  of  beach  sand. 

Simulations  of  a  bed  of  250  identical  spheres  with  radius  R  =  2.5  mm  were  per¬ 
formed  in  exactly  the  same  manner  as  the  composite  particles,  and  the  critical- angle 
value  was  found  to  be  26°  (figure  3),  in  good  agreement  with  experiments  performed 
on  rough  glass  beads  in  rotating  drums  (Jaeger  et  al.  1989). 

(b)  Sheet- flow  simulations 

Both  spherical  and  composite  particles  were  used  in  simulations  of  sheet- flow  trans¬ 
port  to  allow  direct  comparisons  with  the  laboratory  data  of  King  (1991).  Figure  4 
depicts  net  bed-load  transport  rates  from  simulations  with  composite  particles  and 
spherical  particles  along  with  King’s  experimental  data.  Each  point  in  the  plot  repre¬ 
sents  the  average  of  at  least  four  trials  in  both  the  simulations  and  experiments.  For 
the  simulations,  five  consecutive  wave  periods  were  run  with  the  data  from  the  first 
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Figure  4.  The  normalized  flux  is  plotted  against  the  local  bed  slope  for  simulations  with 
spheres  (•),  simulations  with  composite  particles  (x),  and  experimental  data  (o).  All  points 
represent  the  average  flux  from  at  least  four  half-wave  cycles  of  a  sine  wave  with  a  period, 
T  =  4.64  s,  and  maximum  velocity  amplitude  of  1.18  m  s_1.  All  points  are  normalized  by  the 
flux  value  from  the  laboratory  measurements  over  a  horizontal  bed.  The  range  of  values  for 
each  point  typically  does  not  exceed  the  size  of  the  symbol  shown.  The  thin  line  represents  the 
dependence  of  energetics  formulae  for  bed-load  transport  on  the  bed  slope,  tan  /?.  The  thin  line 
is  scaled  to  have  a  value  of  1  for  a  horizontal  bed  in  order  to  isolate  the  dependence  of  the 
transport  formulae  on  bed  slope. 


wave  period  discarded  due  to  transient  start-up  effects.  The  range  of  transport  rates 
for  the  last  four  wave  periods  did  not  exceed  the  size  of  the  symbol  shown.  Likewise, 
for  the  experimental  data  the  range  of  transport  rates  only  slightly  exceeded  the  size 
of  the  symbol  shown  for  two  of  the  points. 

The  time-averaged  transport  rates  from  simulations  with  composite  particles  were 
in  good  agreement  with  the  experimental  data.  Simulations  using  composite  particles 
showed  improved  skill  in  predicting  transport  rates  over  simulations  using  spheres;  in 
all  cases  the  composite  particle  simulations  showed  lower  net  transport  rates  for  the 
range  of  conditions  simulated.  For  the  horizontal  bed  case,  the  average  bed-load  flux 
in  the  composite  particle  simulations  was  within  3%  of  the  average  value  given  by 
the  experiments,  whereas  the  spherical  particle  simulations  were  approximately  25% 
greater  than  the  experimental  value.  For  the  highest  transport  rates  the  composite 
particle  simulations  over-predicted  the  flux  from  the  experiments  by  less  than  25% 
compared  with  60%  over-prediction  for  the  spherical  simulations. 

Although  the  experimental  data  of  King  (1991)  only  provide  bulk  transport  rates, 
the  simulations  allow  detailed  examination  of  sheet-flow  phenomena.  The  average 
time-series  of  the  flux  from  the  four  half-wave  cycles  over  a  horizontal  bed  for  both 
spheres  and  composite  particles  is  shown  in  figure  5.  The  time-series  indicate  that 
the  times  of  initiation  and  cessation  of  motion  are  only  slightly  affected  by  changing 
the  particle  shape.  However,  the  peak  flux  for  the  composite  particle  simulations  is 
approximately  15%  less  than  the  peak  flux  for  the  spherical  simulations.  Integrated 
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over  the  entire  half-wave  cycle  the  net  flux  for  the  composite  particle  simulations  is 
approximately  20%  less  than  the  net  flux  from  the  spherical  simulations. 

The  sheet-flow  simulation  results  were  not  sensitive  to  the  particular  compos¬ 
ite  particle  shape.  We  performed  an  analogous  suite  of  simulations  with  the  ratio 
R2/Ri  =  §,  and  a  single  shape  was  found  to  exhibit  a  critical  angle  similar  to  beach 
sand  by  systematically  varying  the  value  of  72.  When  used  in  the  sheet-flow  sim¬ 
ulations,  nearly  identical  results  were  obtained  to  those  outlined  above,  suggesting 
that  the  simulation  results  are  robust  over  a  range  of  composite  particles  with  shape 
factors  similar  to  beach  sand. 


4.  Discussion 

The  results  from  DEM  simulations  show  how  particle  shape  directly  influences  the 
angle  at  which  a  pile  of  grains  will  fail,  along  with  bulk  bed-load  transport  rates. 
The  simulations  were  conducted  over  a  range  of  conditions  involving  different  par¬ 
ticle  shapes  and  various  bed  slopes.  A  strong  dependency  between  the  grain  shape 
and  bed-load  flux  was  observed.  The  DEM  will  allow  future  work  to  include  more 
simulations  to  elucidate  the  role  of  grain  shape,  particularly  the  degree  of  angularity, 
in  bed- load  transport. 

The  initial  simulations  (figure  3)  demonstrated  that  as  the  particle  shape  was 
changed  from  a  sphere  to  a  dumbbell  the  critical  angle  increased  significantly.  In 
particular,  the  composite  particles  with  R2/R1  =  \  and  72  =  57r/8  have  a  criti¬ 
cal  angle  of  approximately  33°  with  CSF  =  0.84.  These  values  are  characteristic 
of  typical  beach  sand;  this  shape  was  thus  used  in  all  sheet-flow  simulations.  The 
sheet-flow  simulations  using  composite  particles  showed  greater  skill  in  predicting 
the  bulk  bed-load  transport  rates  measured  by  King  (1991)  than  simulations  using 
spherical  particles  (figure  4).  At  small  bed  slopes  and  for  lower  transport  rates  the 
simulations  using  composite  particles  describe  the  relevant  physics  of  the  laboratory 
experiments  using  beach  sand  well.  At  larger  bed  slopes  and  for  higher  transport 
rates  the  simulations  using  composite  particles  showed  better  accuracy  than  spheres 
but  did  not  closely  match  the  experimental  data. 

One  possible  reason  for  the  improved  accuracy  of  the  composite  particles  over 
spheres  is  the  additional  inter-granular  contacts  that  occur  within  a  bed  of  composite 
particles.  For  instance,  any  two  ellipsoids  can  share  only  one  contact  point,  whereas 
the  composite  particles  used  in  our  sheet-flow  simulations  may  share  up  to  four 
contact  points  between  two  adjacent  composite  particles.  In  this  sense,  composite 
particles  offer  a  better  approximation  to  the  shape  of  a  sand  grain  than  other  regular 
shapes  such  as  ellipsoids.  A  bed  of  composite  particles  allow  a  greater  resistance  to 
motion  that  should  reduce  bed  creep  along  with  the  number  of  particles  likely  to 
participate  in  bed-load  transport. 

An  alternative  explanation  is  that  composite  particles  have  greater  rotational  iner¬ 
tia  than  spheres  with  the  same  mass.  The  only  difference  between  simulations  was 
the  particle  shape,  with  the  driving  force  and  governing  equations  for  translational 
motion  identical  between  the  two  sets  of  simulations.  The  peak  flux  (figure  5)  is 
reduced  by  nearly  15%  in  the  time-series  for  the  composite  particles  compared  with 
the  time-series  for  the  spherical  particles.  The  observed  reduction  in  peak  flux  in  time- 
series  for  the  composite  particles  over  spheres  (figure  5)  indicates  that  the  composite 
particles  dissipate  more  translational  kinetic  energy  than  the  spherical  particles. 
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wave  phase 

Figure  5.  The  average  time-series  of  the  bed-load  flux  from  four  successive  half-wave  cycles  in 
simulations  for  spheres  (•)  and  for  composite  particles  (x),  where  R2/R1  —  72  =  57 r/8.  There 

were  no  time-series  of  flux  measurements  from  the  experimental  data  to  be  used  for  compari¬ 
son  here.  The  composite  particle  time-series  showed  deviations  with  initiation  and  cessation  of 
motion  from  the  spherical  time-series.  The  difference  in  peak  flux  between  the  two  time-series 
illustrated  a  bulk  property  of  the  granular  flow  that  the  composite  particles  exhibited  a  larger 
coefficient  of  kinetic  friction  than  the  spherical  particles. 


At  greater  slopes,  the  simulations  with  both  spheres  and  composite  particles  do 
show  significant  deviations  from  the  laboratory  measurements  (figure  4).  The  bed¬ 
load  flux  from  simulations  with  spheres  is  only  in  agreement  with  the  laboratory 
measurements  for  the  smallest  transport  rates.  The  bed-load  flux  from  the  composite 
particle  simulations  remains  in  agreement  with  the  King  experiments  for  the  range 
of  bed  slopes  |tan/?|  <  0.1,  but  over-predicts  flux  for  the  largest  transport  rates. 
Perhaps  the  trend  could  be  explained  by  considering  the  roundness  of  the  particles. 
King  (1991)  reports  that  the  roundness  of  the  quartz  sediments  used  in  the  exper¬ 
iments  were  sub-angular  to  angular,  as  defined  by  Shepard  &  Young  (1961).  Using 
the  same  classification  our  spherical  particles  were  classified  as  well  rounded,  while 
our  composite  particles  were  classified  as  rounded.  Composite  particles  have  greater 
rotational  inertia  than  spherical  particles  with  the  same  mass.  Likewise,  angular 
grains  have  greater  rotational  inertia  than  composite  particles  with  the  same  mass. 
The  rotational  inertia  of  individual  particles  as  their  shape  goes  from  having  three 
axes  of  symmetry  (spheres)  to  one  axis  of  symmetry  (composite  particles)  to  no 
symmetry  axes  (angular  grains)  will  increase  bulk-energy  dissipation  with  increasing 
energy  and  should  result  in  lower  bed-load  flux  for  the  most  angular  sediment  grains. 
The  results  of  the  simulations  suggest  that  particle  shape  is  an  important  parameter 
for  predicting  bed-load  flux. 
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(a)  Sediment  shape  dependence  in  bed-load  transport  formulae 

Energetics  models  (Bowen  1980;  Bailard  &  Inman  1981)  describing  bed-load  trans¬ 
port  in  the  surf  zone  based  on  BagnokTs  (1966)  description  of  bed-load  transport  in 
unidirectional  flows  contain  only  one  parameter,  tan(</>),  that  depends  on  the  shape  of 
the  sediment  particles,  where  <j)  is  the  internal  angle  of  static  friction  of  the  sediment. 

We  have  shown  that  a  composite  particle  can  be  constructed  to  exhibit  a  value  of  <j> 
similar  to  beach  sand.  Now,  we  can  examine  more  closely  the  particle  shape  depen¬ 
dence  in  energetics  models  for  bed-load  transport.  For  the  simulated  conditions  of  a 
unidirectional  half-period  sine  wave  the  bed- load  formula  of  Bailard  &  Inman  (1981) 
depends  on  tan  <fi  as  follows: 


where  q  is  the  bed-load  transport  rate  and  tan/3  represents  the  local  bed  slope  along 
the  flow  direction. 

For  flow  over  a  horizontal  bed,  tan  /3  =  0,  the  bed-load  transport  rate  is  inversely 
proportional  to  the  coefficient  of  static  friction,  tan  </>.  The  time-averaged  bed-load 
flux  over  a  horizontal  bed  for  simulations  with  composite  particles  is  in  excellent 
agreement  with  the  experimental  data  for  the  single  point  shown  in  figure  4.  The 
simulations  with  spheres  gave  an  increase  in  bed-load  flux  of  approximately  25%  over 
the  King  experiments  for  the  case  of  the  horizontal  bed.  The  value  of  tan  <j>  for  the 
spheres  in  simulations  was  found  to  be  about  0.5.  The  value  of  q  in  equation  (4.1)  is 
approximately  25%  greater  when  tan</>  —  0.5  than  when  tan</>  =  0.63,  with  tan/3  =  0, 
suggesting  that  for  bed-load  transport  over  a  horizontal  bed,  the  shape  parametriza- 
tion  provided  by  energetics  models  is  consistent  with  both  spherical  and  composite 
particle  simulations. 

For  sloping  beds,  the  bed-load  transport  rate  given  by  equation  (4.1)  is  a  linear 
function  of  the  local  bed  slope,  tan/3.  The  thin  line  shown  in  figure  4  represents  the 
dependence  of  bed-load  transport  on  local  bed  slope  given  by  equation  (4.1)  with 
tan</>  =  0.63  (Bagnold  1956).  Since  King  (1991)  did  not  report  the  value  of  tan</> 
for  the  grains  used  in  his  laboratory  experiments,  the  vertical  offset  of  the  line  was 
adjusted  so  that  the  flux  value  for  a  horizontal  bed  goes  through  the  value  of  the  flux 
from  the  laboratory  data.  The  laboratory  data  are  in  good  agreement  with  the  ener¬ 
getics  dependence  for  the  entire  range  of  bed  slopes  shown,  given  the  choice  of  tan  <f>. 
The  composite  particle  simulations  are  also  in  good  agreement  for  the  range  of  bed 
slopes  —0.1  <  tan/3  <  0.1.  Typically,  it  is  not  common  to  find  bed  slopes  in  the  surf 
zone  much  greater  than  about  0.1  except  on  the  slip  faces  of  bed  forms,  suggesting 
that  our  composite  particle  model  captured  the  essential  physics  for  particle  shape 
and  local  bed  slope  parametrized  by  energetics  models  in  equation  (4.1).  Our  expec¬ 
tation  is  that  for  higher  slopes  where  flux  values  can  be  large,  particle  shape  will 
have  a  substantial  influence  on  bed-load  flux  mechanics,  necessitating  a  modification 
to  the  formula  given  in  equation  (4.1). 

5.  Conclusions 

We  have  presented  a  new  DEM  for  sheet-flow  transport  of  coarse  sediments  in  a 
WBBL  under  oscillatory  flow.  The  DEM  used  composite  particles  to  simulate  the 
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non-spherical  shape  of  beach  sand.  Simulations  performed  to  measure  the  internal 
angle  of  static  friction  for  beds  of  composite  particles  reveal  that  there  exists  a  com¬ 
posite  particle  that  exhibits  a  similar  coefficient  of  static  friction  to  beach  sand, 
tan<£  «  0.63.  Additionally,  this  composite  particle  has  a  CSF  slightly  greater  than 
about  0.8,  in  the  upper  range  for  CSF  values  of  naturally  occurring  sediments.  This 
composite  particle  shape  that  most  closely  matches  the  properties  of  beach  sand  was 
then  used  to  simulate  sheet-flow  transport  conditions  observed  in  the  laboratory  data 
of  King  (1991).  The  simulation  results  showed  that  the  bed-load  flux  for  gently  slop¬ 
ing  beds,  |tan/?|  ^  0.1,  with  sinusoidal  wave  forcing,  agrees  well  with  the  laboratory 
measurements.  However,  the  simulations  suggested  that  particle  roundness  could  be 
an  important  parameter  for  predicting  bed-load  flux  at  larger  bed  slopes  and  higher 
transport  rates. 

The  simulation  results  suggest  that  the  composite  particle  approximation  to  the 
shape  of  a  sand  grain  is  useful  and  robust  for  discrete-element  modelling  of  bed-load 
transport  phenomena.  Comparisons  of  simulations  with  both  spheres  and  composite 
particles  with  laboratory  data  indicate  that  particle  shape  is  an  important  variable 
for  understanding  the  bulk  properties  of  sheet  flow. 

J.C.  thanks  P.  K.  Haff  for  valuable  discussion  motivating  this  research.  We  thank  C.  S.  Thaxton 
for  a  review  of  the  analytical  solution  provided  in  the  paper.  This  work  was  performed  while 
J.C.  held  a  National  Research  Council  Research  Associateship  Award  at  the  Naval  Research 
Laboratory.  Research  was  supported  by  the  Office  of  Naval  Research. 
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